First-passage and first-exit times of a Bessel-like stochastic process 
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We study a stochastic process Xt related to the Bessel and the Rayleigh processes, with various 
applications in physics, chemistry, biology, economics, finance and other fields. The stochastic 
differential equation is dX t — (nD/X t )dt + \/2DdWt, where Wt is the Wiener process. Due to the 
singularity of the drift term for X t = 0, different natures of boundary at the origin arise depending 
on the real parameter n: entrance, exit, and regular. For each of them we calculate analytically and 
numerically the probability density functions of first-passage times or first-exit times. Nontrivial 
behaviour is observed in the case of a regular boundary. 
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I. INTRODUCTION 



In the theory of stochastic processes, the first-hitting 
time is defined as the time when a certain condition is ful- 
filled by the random variable of interest for the first time; 
it is random itself and a particular case of a stopping 
time. We speak of a first-passage time when the random 
variable reaches a certain level for the first time, and of a 
first-exit time when it leaves a certain interval for the first 
time. A standard example of a first-passage time prob- 
lem is the decision of an investor to buy or sell a stock 
when its fluctuating price reaches a certain threshold. 
However, first-passage times play an important role also 
in chemical physics; early examples are given by mod- 
els describing the dissociation of diatomic molecules as 
a first-passage time problem, where dissociation occurs 
if a certain critical energy level is reached through colli- 
sions |l|-|3J]. A view on diffusion in fluids based on first- 
passage times has been proposed by Munakata [4j] , where 
self-diffusion is measured via the first-passage time with 
respect to a boundary marked by a sphere centered at 
the original position of a labeled particle. Problems like 
neuron dynamics, self-organized criticality or dynamics 
of spin systems can be viewed as first-passage processes 
in one dimension [5| . The first-passage problem is closely 
connected to persistence, which is the probability that a 
random variable does not leave a certain region up to 
a certain time, i.e. the complementary event to a first- 
passage at the same time. The problem of persistence in 
spatially extended nonequilibrium systems has attracted 
great interest both theoretically and experimentally, see 
Majumdar [6| and references included therein, where per- 
sistence is defined as the probability that for an arbitrary 
nonequilibrium field 0(r, t) the quantity <p(r, t) — (<p(r, t)} 



does not change sign. The nonequilibrium field can also 
be a scalar or tensorial order parameter field. Yurke et 
al. [3] have measured the probability that the local or- 
der parameter in a twisted nematic liquid crystal sys- 
tem has not switched its state up to a time t. Persis- 
tence phenomena have also been studied in the context 
of phase-ordering dynamics @, M, diffusion fields [1C| . 
and reaction-diffusion systems [fl(. All these systems 
share the characteristic property that persistence, and 
hence also the distribution of first-passage times, follows 
a power law with some non-trivial exponent. However, in 
the literature that we have read, the reference point with 
respect to which persistence was measured has always 
been zero. In this work we shall consider the first-passage 
or first-exit problem with respect to a certain level b for 
a stochastic process that may or may not be able to cross 
the origin, depending on the nature of the boundary at 
zero. 

Our model can be described by the stochastic differen- 
tial equation 



dX t = ^dt 

Xt 

where W t is the Wiener process with zero mean 

(W t ) = 0, 
and the autocovariance function 

(W t Wt>) =xmn(t,t'), 



(1) 



(2) 



(3) 
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the constant diffusion coefficient D is positive, and the 
real parameter n controls the relative strength of the drift 
and diffusion terms of the model. 

Except for n — the origin X t = is a singular point, 
the nature of which depends on the value of n. Intuitively 
one might think that the stochastic process cannot cross 
the origin for a non-zero n: it should be bounded to the 
interval (0, oo) or (— oo, 0) depending on the initial value 
xq of the process. As we shall see later, this is only true 
for a certain range of n. 



For n — the process is nothing but the Wiener pro- 
cess. The probability density function (PDF) of the first- 
passage time of a certain level b at time T starting at xq 
is well known, 



f(T) = 



x 



T-3/2 



V4vr£> 



■ exp 



4DT 



(4) 



This result can be derived, for example, in an elegant 
way by a simple scaling argument pj]. For long times 
one obtains from Eq. ((4]) a power law f(T) ex T~ 3 / 2 . 

For nonzero n the situation is not so simple any more. 
For n < 1 Bray [13| obtained a result for the PDF of 
the first time to hit the origin, which for long times 
is a power law too, f(T) ex T~( 3- ")/ 2 ; see also below. 
The persistence probability is then simply the probabil- 
ity that no hit has occurred in the time span T, that is 
1 - f T /(T') dT' = /" f(T') dT' and is again a power 
law, with exponent —(1 — n)/2. 

In this paper we consider a more general problem. We 
ask for the first time to leave the interval (0, b) either by 
crossing the upper boundary b or by hitting the origin, if 
the nature of the singularity allows the latter event. Cor- 
respondingly we calculate the first-passage time distribu- 
tion with respect to the upper bondary b or the first-exit 
time distribution for the interval (0,6). 

The PDF is determined as the solution of a bound- 
ary value problem, which is a Sturm-Liouville eigenvalue 
problem. On the semi-infinite interval (0, oo) the problem 
has a continuous spectrum, whereas on the finite interval 
(0, b) the spectrum is discrete. 

The paper is organized as follows. In Sec.|TT]we shortly 
explain the scientific relevance of our process and discuss 
its relation to other model processes. The main results 
of the paper are presented in Sees. IIIII and IIV1 In Sec. Mil 
we analyze the nature of the singular point at the ori- 
gin, which depends on the value of n. First we adopt 
a heuristic approach by Bray [l3[, then we present a 
more sophisticated analysis following a scheme proposed 
by Feller [3|- In Sec. [TV] we calculate the PDFs of the 
first-passage time or of the first-exit time; we derive the 
backward Kolmogorov equation, formulate the boundary 
value problem, which is of Sturm-Liouville type, and give 
the general solutions for different ranges of n. We con- 
clude Sec. IIVI giving a short description of a numerical 
simulation method and comparing the analytical results 
with those from simulation. 



II. PHYSICAL MOTIVATION AND RELATED 
MODELS 

First we remark that Eq. ([T]), setting n = d— 1 — U/D, 
governs the dynamics of the radial component of the posi- 
tion of a random walker in a logarithmic potential U log x 
in dimension d or of a free random walker in an effective 
dimension d' = d — U / D, which may have noninteger 
values |13| . In this context it appears natural to assume 
x > 0, and for free diffusion U = to restrict to n > —1. 



Eq. ([T]) appears in various physical, chemical and bio- 
logical problems. The context in which the relevance of 
Eq. (H|) arises will now be explained starting from generic 
considerations and then proceeding further_with specific 
physical problems. Godreche and Luck [15| introduced 
a classification of stochastic processes into a group with 
"narrow" distributions, where all moments are finite, and 
"broad" distributions, where PDFs exhibit a power-law 
decay, and hence only a finite number of moments con- 
verge. The PDFs of the persistence, and therewith also 
of the first-passage time, will decay respectively either 
faster than any power law or algebraically, depending on 
the nature of the process imposed by its distribution. 

The motion of atoms in a one-dimensional optical lat- 
tice formed by two counterpropagating laser beams with 
linear perpendicular polarization was studied by a similar 
equation in the high momentum region, where the mo- 
mentum takes the role of the stochastic variable x (16l | . 

The Barkhausen noise was described phenomenologi- 
cally by a model where the domain wall velocity as a 
function of the magnetization is also described by a sim- 
ilar Langevin equation if the demagnetizing factor is ne- 
glected [13, LL8( • The magnetization takes the role of time. 

Fogedby and Metzler, the former of which had alr eady 
analyzed the generic Langevin equation before [13, [20] , 
have applied the model to study the variable size of a 
DNA "bubble" [2l|, which emerges when at a certain 
temperature hydrogen bonds connecting base pairs from 
the opposite strands of the double helix are broken. The 
bubble size is measured by the number of broken bonds; 
in a continuum limit the discrete number of broken bonds 
can be replaced by a continuous variable x and, according 
to the Poland-Scheraga model [22], the free energy of the 
system can be approximated for small bubble sizes as 



T sa cfceTloga;, 



(5) 



where c is a positive constant. Equilibrium is reached for 
a minimum of the free energy and the dynamics follows 
the Langevin equation 



dx 
~dt 



-D 



dT 
dx 



m, 



(6) 



where t;(t) is a thermal noise assumed to be Gaussian 
with autocovariance (£(*)£(*')) = 2Dk B TS(t - t'). 

A similar model was employed studying the translo- 
cation of a polymer through a pore, where the number 
of monomers on one side is chosen as the "translocation 
coordinate" (23J . 

An interesting application of the model was found by 
Bray [131 ] . who showed that persistence and nonequi- 
librium critical dynamics are related in the context of 
the two-dimensional AY"-model with non-conserved or- 
der parameter, where the critical temperature is the tem- 
perature Tkt of the Kosterlitz-Thouless phase transition. 
The dynamics of a vortex-antivortex pair can be mapped 
to a one-dimensional Langevin equation corresponding to 
Eq. (fTJ) by a series of transformations. 



It is impossible for us to discuss comprehensively all 
the publications dealing with similar models because too 
many of them exist. The list becomes even longer con- 
sidering Eq. ([T]) as a special case of several different more 
general types of stochastic process. 

On the one hand it is a special Rayleigh process [24|, 



dX t = (n- X Xt l +ii x X t )dt- 



rdW t 



(7) 



where /U_i, [i\ and a are constants. Setting /i_i = nD, 
fj,i = and a — \J2D results in Eq. ([1}, while setting 
/i_i = reproduces the radial Ornstein-Uhlenbeck pro- 
cess in one dimension [20] . Using methods of classical Lie 
group symmetry analysis, it was shown that the Rayleigh 
process belongs to a maximal invariance group with six 
parameters whose member equations can be reduced to 
the standard diffusion equation by appropriate changes 
of variables, thus leading to the analytical solutions of 
these equations [27H29J . 

The Rayleigh process is widely employed in physics, 
economics, finance, and other fields, e.g. biometry [30|. 
In physics, Eq. (J7|) appears e.g. within non-abelian gauge 
theories in the framework of stochastic quantization [3l| . 
In economics, Eq. ([7]) appears e.g. as a special case of a 
more general diffusion process whose stationary solution 
has been proposed to model the distribution of the profit 
rate of firms [32j. In finance, the applications usually 
employ the form of the generalized Bessel process, which 
is introduced below. 

On the other hand our process can be mapped onto 
the Bessel process 



dY t = adt + b-s/YtdWt 



(8) 



via the transformation Y t — X 2 : multiplying Eq. (TTJ) 
with X t , interpreting the stochastic integral in the Ito 
sense and using Ito's lemma, one recovers Eq. (|5|) with 
a = 2^_i +<j 2 = 2(n + l)D and b = 2a = V^D, where 
we have assumed Y t > 0. 

In finance, extensions of the Black-Scholes-Merton 
(BSM) option pricing formula [33, [3J] based on diffusion 
processes where the volatility is a function of the under- 
lying, called constant elasticity of variance or Cox pro- 
cesses [3a], reduce to a more general Bessel process with 
an additional term proportional to Y t in the drift, corre- 
sponding to the Raleigh process, by means of a non-linear 
transformation and a measure change. The generalized 
Bessel process 



dY t = (a + aiY t )dt- 



bYfdWt 



(9) 



describes the underlying stock price in the BSM lognor- 
mal model with /3 = 1, ao = and a\ > 0, the short 
interest rate in the Vasicek model [36| with j3 = 0, ao > 
and a\ < 0, and the short interest rate in the Cox- 
Ingersoll-Ross short rate model (37| with /3 = 1/2, ao > 
and a\ < (when a\ < the process is called mean- 
reverting). These models with three adjustable parame- 
ters are all solvable. A general solution was proposed for 



a larger family of models with up to seven parameters; 
it has a similar structure as the BSM formula, the most 
notable difference being that error functions are replaced 
by confluent hypergeometric functions [381 ] . 

Eq. ((9]) is used to describe the underlying with (3 = I 
or f3 — 1/2 and ao — also when pricing path-dependent 
options, e.g. barrier and lookback 39-42] or Asiatic op- 



tions [43j. In this context the cumulative probability 
distribution F(T) of the first-passage times of an upper 
barrier, i.e. the probability that the barrier is reached 
within a time T, is the probability that an up-and-in, or 
knock-in, option is valid at its maturity T (here the bar- 
rier is an entry point), while I — F{T) is the probability 
that an up-and-out, or knock-out, option is valid at T 
(here the barrier is an exit point). In both cases a valid 
barrier option behaves as a European option; thus F(T) 
is the probability that at maturity an up-and-in option 
behaves as a European option, and 1 — F(T) is the prob- 
ability that at maturity an up-and-out option behaves 
as a European option. These considerations lead to one 
approach (among others) to price barrier options. 



III. THE NATURE OF THE SINGULAR POINT 
AT THE ORIGIN 

The quantity we are interested in is the first-passage 
time with respect to a certain level b, when the initial 
value x satisfies < x < b. Clearly, the upper limit of 
the process is the artificially set absorbing boundary at 
x = b. It requires some effort to understand the nature 
of the lower boundary x = 0, which is a singular point of 
the stochastic differential equation. We shall first adopt 
the heuristic argumentation by Bray [f3j and then apply 
a more sophisticated classification scheme proposed by 
Feller El. 



A. Heuristic arguments 

The Fokker-Planck equation corresponding to the 
stochastic differential equation (fTJ) is 



dpjx, t) 
dt 



d_ 
dx 



nD 

P{x,t) 

X 



D 



d 2 p(x, t) 
dx 2 



(10) 



This does not depend on whether Eq. (fTJ) is interpreted 
in the Ito or Stratonovich sense because in our case the 
noise is additive, i.e. independent of x J44(. We restrict 
to s > 0, the case x < being symmetric to the pre- 
vious one. The general solution of the Fokker-Planck 
equation (|10p requires the knowledge of two linearly in- 
dependent solutions. Using the separation ansatz [l3j 



p{x,t)=x < - 1+n ^ 2 R k (x)e- 



Dk z t 



one gets the Bessel differential equation 



d 2 R k , ldR k 



dx 2 



x dx 



k 2 -—)R k = U. 



(11) 



(12) 



where v — (1 — n)/2, whose solutions are the Bessel 
functions of the first kind J„ (kx) and of the second kind 
Y v {kx). 

For non-integer v also J u {kx) and J- V {kx) are linearly 
independent, and we can use J- v instead of Y v . In this 
case the general solution of Eq. (fT0|) can be written as 



p{x, t) = x 1 -" / [A(k) J v {kx) + B{k)J_ v (kx)} 
Jo 

•xe- Dk ' 2t dk, (13) 

where the coefficients A(k) and B(k) are to be deter- 
mined by initial and boundary conditions. 

The Bessel function of the first kind is given by [45| 



J v {kx) = 2J 



(-1)' 



/=o 



k 1 



UT(l + v + l) V 2 



21+v 



(14) 



Denoting the contributions to p(x, t) coming from J v and 
J- v by pu(x,t) and p_„(a;, £) respectively, we can use 
Eq. (Q3]) to write Eq. (jTHJ in the form 



p(a;,t) = p v {x,t) + p^ v (x,t) 

00 

= Y,[cl(t)x 2 ^+cl u (t)x 21 ^], (15) 



;=o 



with 



, l\T(l + u + l) 



c l u (t) = / A(fc) 



(-l)<(fc/2) 2 ^ n fe2i 



c -^=/ B <*hm-v + i) 



dk 

< - "dk. (16) 



The behaviour of the PDF for x — >• + is determined by 
the leading order terms and thus, approaching zero, we 
find 



p(x,t) ~ cl{t)x + c°_ u (t)x Tl 



(17) 



For 1 < n the leading order term is cP v {t)x and the next 
to leading order term is <P_ v (t)x n . For —1 < n < 1 it 
is vice versa. For n < — 1 the contribution from J_„ 
has a non-normalizable singularity at the origin; see also 
below. Introducing the probability current density 



n d 
j(x,t) = D( -- — )p( x ,t), 



(18) 



the Fokker-Planck equation (|10|) can be written in the 
form 



dtp + d x j = 0. 



(19) 



Terms proportional to x n do not contribute to j(x,t), 
and so we arrive at 

j{x,t)~<Z(t)D(n-l), (20) 

which holds in leading order for x —¥ + . 



For 1 < n the coefficient c°(i) is positive because in 
this case c°(t) = lxm x _ > .Q+p(x,t)/x, and correspondingly 
j(x,t) > near the origin. 

For — 1 < n < 1 this is not true in general, see Eq. (jTTJ). 
The characterization can be made clearer imposing spe- 
cific boundary conditions. 

For example, absorbing boundary conditions at x = 
require lim 2 ,_ > o+ p{ x i t) = an d H 1 oi x ^q+ j(x,t) < 0. 
For — 1 < n < the former condition can be fulfilled 
only setting B(k) — 0. A(k) is determined by the initial 
condition. Observing the orthogonality relation [46| 

6(a-0) = a U v {ak)J v {0k)dk, (21) 

Jo 

which holds for v > —1/2, i.e. for n < 2, we see that 
choosing A(k) such that 



p(x,t)=x 1 -"x% kJ v {kx )J v {kx)e- Dk * dk (22) 

Jo 

fulfills the initial condition p(x, 0) = S(x — Xo)- 

The integral in Eq. (|2"21 can be explicitly evaluated [47| 
with the result 



P(,.0=,'-^4«p(-4±^)7„(=l), 



(23) 



given already by Bray [13| . Here I v is the modified Bessel 
function of the first kind defined by 



w = E 



1 



k=0 



k\T(k + u + l) 



n(f) 



£ \ i/+2fe 



(24) 



Since in this case -B(a;) = 0, which implies that also 
c_„ = 0, we have \ivci x ^ + p(x,t)/x = c v (t) > 0, and 
from Eq. (|2"0)l it follows that j(0 + ,t) < as required. 

Note that the case of free diffusion (n = 0, i.e., v = 
1/2) with an absorbing boundary condition imposed at 
the origin is included in Eq. (|2"3")l . A short calculation 
gives the well known result which can be obtained, e.g., 
by the mirror method. 

As already shown by Karlin and Taylor [48| , for n < 1 
total absorption at the origin occurs in finite time. Cor- 
respondingly, in this case there exists a stationary solu- 
tion of the Fokker-Planck equation which is a Dirac delta 
function S(x); see also Alfarano et al. [32(. Formally this 
can be seen as follows. Observe that Eq. (fj"0|) admits sta- 
tionary solutions p s (x) ex x n whith n < 1 that are not 
normalizable at the origin. Using the concept of weak 
normalization introduced by Senf et al. 49], it can be 
shown that the weakly normalized version of p s {x) is just 
a Dirac delta distribution, pj(x) — S(x), in the sense 
that LpJ (x)ip(x)dx — tp(xo), where ip(x) is a test func- 
tion and xq is included in the support S. 

More insight into the qualitative behaviour of the sys- 
tem near and at the origin is provided by the classification 
scheme of Feller which is discussed in the next subsection. 



B. Formal classification 

The modern classification of the boundaries of diffusion 
processes has been developed by Feller [f4| and is based 
on semigroup operator arguments. We shall now briefly 
review the necessary theory for the boundary classifica- 
tion employing the notation of Karlin and Taylor [48| in 
order to be able to classify the origin for our process. 

In the following let X t be a process defined on the 
interval / = (I, r), where the two endpoints can be both 
finite or infinite. Also let the process start at the initial 
value Xq = x, and a and b be two finite real numbers 
such that the inequality l<a<x<b<r holds. We 
shall consider regular diffusion processes in the interior of 
/, i.e. processes for which the first-passage time T y with 
respect to an arbitrary level y in the interior of / is finite 
with a positive probability 



P(T y <oo|A = x) > 0. 

The three central quantities are 

u{x) = P{T b < T a \X = x), 
v(x) = (T*\X Q = x), 

w ( x ) = ( / g(X s )ds\X = x 



(25) 



(26) 

(27) 

(28) 



where g is an arbitrary functional of the stochastic pro- 
cess, and we have defined T* = T a ^ = mm{T a ,Tt,}. It 
can be shown |48] that under certain conditions these 
quantities satisfy the boundary value problems 

Lu(x) = 0, u(a) = 0, u(b) = I, (29) 

Lv(x) = -1, v(a) = 0, v(b) = 0, (30) 

Lw(x) = -g(x), w{a) = 0, w(b)= 0, (31) 

with the differential operator L acting on a function f(x) 
as follows: 

Lf(x)=n(x)f(x) + \a 2 (x)f"(x). (32) 

The proof for u(x) invokes the law of total probability 
and uses a Taylor expansion to the second order around 
x of the functional u(Xh) at a small instant of time h. 
The proof for w(x) uses a similar procedure, and finally 
the case v(x) follows as a special case of w(x) by setting 
g(x) = l. 

The differential operator given by Eq. (|32[) can be writ- 
ten as 



Lf{x) = \a\x) ( 2 4f\f'(-) 
2 \er z (a;) 



n*: 



i 



2m(x) dx 



1 df(x) 



s(x) dx 



with 



s(x) = cxp 



MO 



ds 



(33) 



(34) 



(the lower integration boundary is not indicated because 
it is arbitrary) and the speed density 



m(x) 



1 



a 2 (x)s(x) 
Introducing the scale function 



S(x) 

and the speed function 

M(x) 



(35) 



(36) 



(37) 



(38) 



The definitions given by Eqs. (I34H37[) naturally induce 
measures of closed intervals J = [c, d] : the scale measure 



s(r)) di] 



m(r)) dr/, 



Eq. (|33[) can be rewritten in the form 



Lf{x) = - 



I d 



2 dM{x) 



df{x) 



dS{x) 



S[J] = S[c, d] = S(d) - S(c) = J s(x) dx, 



(39) 



and the speed measure 



M[J] = M[c, d] = M(d) - M(c) = / m(x) dx. (40) 

J c 

These measures are fundamental for the classification of 
diffusion processes. The scale measure for an infinites- 
imal interval J = [x, x + dx] is written symbolically as 
S[dx) = S(x + dx) — S(x) — dS(x) — s(x)dx, and the 
same applies for the speed measure. 

Then Eqs. ([29]) and (|3Tj) can be easily integrated first 
with respect to the speed measure and thereafter with 
respect to the scale measure. Using the notation intro- 
duced above, the solutions can be expressed in compact 
form as 



u(x) 



S[a 



S[a, b] ' 



(41) 



and herewith 
w(x)=2lu(x)J S[ri,b]g(ri)dM(T,) 

+[!-«(*)] / S[a,r,]g(ri)dM(v) 



(42) 



The solution of Eq. (|30[) follows again from the special 
case of g(x) = 1 in Eq. (|42|) . 

In the following only those definitions relevant for our 
classification will be mentioned, and not every proof can 
be given in detail. The book by Karlin and Taylor [48[ is 
excellent for a deeper understanding. For the classifica- 
tion of the left boundary / of a process, the procedure is 
to regard u(x) and v(x) in the limit a —¥ I. An analogous 



approach is employed for the right boundary; however 
we shall only be interested in the left boundary, which in 
our case is the zero level. 

The first definition which is important for the under- 
standing of whether a boundary can be reached is the 
attractiveness. A left boundary is called attractive if 
S(l,x) := lim _>/ S[a, x] < oo for some x <E (l,r). If 
the scale measure S(l,x] is finite for some x € (l,r), this 
is also true for all x in this interval. Hence it follows di- 
rectly from Eq. flSJ) that P (T/ < T b \X = x) > for all 
I < x < b < r, i.e. there is a positive probability that the 
left boundary is reached before the level b in the interior 
of the interval, provided that the former is finite. 

The next question is whether a boundary is attain- 
able in finite time. This can be measured by lim a _,.; v(x), 
which is the expectation value of the first exit time from 
the interval (l,b). Provided that the boundary is attrac- 
tive, and using the solution v(x) given by setting g(x) = 1 
in Eq. (j4"!2|) , it can be shown that it suffices to check 
whether a certain functional called E(Z) is finite in order 
to establish the attainability of the boundary. Hence a 
left boundary is said to be attainable if it is attractive 
and the functional 



E(0 



J S(l,Z]dM(£) = J M[ v ,x]dS( v ) 



(43) 



is finite, otherwise it is said to be unattainable. Similarily 
one can define 



M(U]dS(0= / S[ri,x]dM(ri). (44) 



N(l) 



The classification of the left boundary of a process is 
based on whether the functionals S(l,x], M(l,x], £(/), 
and N(l) are finite or not. These functionals are not 
independent of each other and some combinations are 
impossible; for example an attainable boundary is always 
attractive. 

Using Feller's terminology, four types of boundaries 
can be distinguished. A process can both enter or leave 
from a regular boundary. The criteria for a left boundary 
to be regular are S(l,x] < oo and M(l,x] < oo. In the 
case of an exit boundary it is impossible to reach any inte- 
rior state b if the starting point approaches I. A boundary 
is an exit boundary if E(Z) < oo and M(l,x] — oo. An 
entrance boundary cannot be reached from the interior 
of the state space, but it is possible to consider processes 
beginning there. It suffices to show that S(l,x] = oo 
while N(l) < oo to prove that I is an entrance bound- 
ary. Finally, a natural or Feller boundary can neither be 
reached in finite mean time nor be the starting point of 
a process, and the corresponding criteria are E(Z) = oo 
and N(l) = oo. 

We are now able to classify the zero level of our pro- 
cess. The first step is to check the attr activity. The pa- 
rameters determining our process are //(x) = nD/x and 
cr 2 (x) = 2D. Since the scaling function only depends on 
the upper integration limit, we can choose the lower limit 



dz I = r/ 



in a convenient way such that 

s(n) = cxp 
Then the scale measure of interest is 
5(0, x] = lim / r]- n di] 

_ J T^ ( xl ~ U ~ lim ^o a 1 "") for n ? 1, 



(45) 



log x - lim a ^o log a 



for n = 1, 



(46) 



and thus the origin is attractive (5(0, x] < oo) for n < 1, 
and non-attractive (5(0, x] = oo) for n > 1. 
The speed density of the process is 



mfa) = ^, 



(47) 



and we can evaluate the speed measure of an interval 
(0, x] as 



M(0,x] = — lim / r] n df] 

2D a— s-0 



2D(n+l) 



„n+l 



lim a _j. a 



n+:n 



^ (log x - lim a ^ log a) 



2D 



for n^-1. 
for n = — 1. 

(48) 



Hence we have M(0, x] < oo for n > —1 and M(0, x] = oo 
for n < — 1. We now have established the nature of the 
zero level for n<l:ifrt<~ 1 the origin is an exit bound- 
ary and in the case — l<n<litisa regular boundary. 
A regular boundary in the origin is the most complicated 
case. Karlin and Taylor [48] describe a regular boundary 
as follows: 

For a regular boundary a variety of boundary 
behaviour can be prescribed in a consistent 
way, including the contingencies of complete 
absorption or reflecting, elastic or sticky bar- 
rier phenomena, and even the possibility of 
the particle (path), when attaining the bound- 
ary point, waiting there for an exponentially 
distributed duration followed by a jump into 
the interior of the state space according to a 
specified probability distribution function. In 
the latter event, the process only exhibits con- 
tinuous sample paths over the interior of the 
state space. 

The last step is to compute N(0) for the classification 
of the case n > 1. Using Eq. (|44l) we get 



N(0) 



S[r),x]dM(rj) 



i) 



*(fc"<)n"<i>i- 



l 

2D 



(49) 



It is easy to show that this double integral is always finite, 
and thus the origin is an entrance boundary for n > 1. 

Summarizing, the nature of the boundary at zero has 
the following behaviour: exit if n G (—00, — 1], regular if 
n <G (— 1, 1), and entrance if n € [1, 00). 



IV. FIRST-PASSAGE AND FIRST-EXIT TIMES 

A. Heuristic approach to first-passage times 

In a somewhat heuristic approach the first-passage (or 
first-exit) time PDF to leave an interval (I, r) can be writ- 
ten as f(T) = -G(T) ;50], where 



G(T) 



p(x, T) dx 



(50) 



is the probability that the particle is at time T in (a, b) 
when it has started at zero time at xq and we have calcu- 
lated the PDF p(x,t) imposing absorbing boundary con- 
ditions on those boundaries where the particle can leave 
the interval. Writing the Fokker-Planck, or forward Kol- 



mogorov, equation in the form d t p 
readily 

G(T) = -j(x,T) 



d x j = we have 



(51) 



If the upper boundary is a natural boundary at 00, 
j(oo,T) — 0, and we are interested in hits at the origin 
we have f{T) — —j(0,T). For our problem the probabil- 
ity current density at the origin has been calculated in 
Sec. lIHAI From Eqs. (J2UJ) and fl23) one obtains 



/CO 



r(i/) \4D 



T 2 



T 



-0+1) 



exp 



ADT 



(52) 



For long times this is a power law f(T) oc T^ 3 -™)/ 2 [H[. 

The result was obtained solving the Fokker-Planck 
equation (|10[) on the semi-infinite interval (0, 00) with 
an absorbing boundary condition at x — and the initial 
condition at x = xq, and it is restricted to n < 1, i.e. 
v > 0. The spectrum of this boundary value problem is 
continuous. 

However, if we are interested in the first time to leave 
a finite interval, we have to solve a boundary value prob- 
lem with, for example, absorbing boundary conditions at 
both ends of the interval which typically has a discrete 
spectrum. We find it preferable to adopt a more formal 
approach, based on the backward Kolmogorov equation. 
The boundary value problem can then be transformed to 
a canonical Sturm-Liouville problem and systematically 
solved. 



B. The backward Kolmogorov equation 

In this section we use a special Fokker-Planck tech- 
nique proposed by Kearney and Majumdar [51j to ob- 
tain a differential equation for the first-passage time PDF 



in Laplace space. Their method is very powerful, be- 
cause the boundary conditions can be easily established 
in Laplace space and the functional V[-Xt] can be chosen 
such that different relevant quantities can be computed. 
Therefore we present the application of this method to 
our problem in some detail. 

Considering a stochastic process starting at Xq = x 
governed by the stochastic differential equation ([1}, we 
are interested in the PDF f(Tf,,x) of the first-passage 
time Tf, with respect to a certain level b, i.e. the time 
when the process has reached the level b for the first 
time. First of all we define an arbitrary functional T^[A^ t ] 
by 



T = 



T h 



V[X t ] dt. 



(53) 



T can have several meanings; in the special case V[Xt] = 
1 it is simply the first-passage time Tf,. The strategy is to 
find a differential equation in Laplace space for f(T,x). 
The Laplace transform of f(T, x) with respect to T is 
given by 

f(s,x) - C T [f(T,x)](s) 

f(T,x)e' sT dT=(e- sT ) T , (54) 



where s G C. Splitting the interval [0,7},] into a small 
interval [0, At] and an interval (At, Tf,], we can expand 
the integral over the small interval to first order in At: 

/ V[X t ]dt = V[x]At + o(At). (55) 

Jo 

Thus Eq. (53|) becomes 



T = V[a:]At+/ V[X t ]dt=:T 1 +T 2 . (56) 

J At 



Inserting Eq. (|56|) into Eq. (|54|) gives 

= / f(T,x)e' sTl e' sT2 dT. (57) 

Jo 

If we split the interval [0, 7&] as described above, we must 
take into account that we also split our trajectory in two, 
where the starting point of the second part, y := X/± t — 
x + Ax, is random itself. Therefore the PDF takes the 
form 



f(T,x) 



f(T 1 ,x)f(T 2 ,y)dy 



b—x 



f(T 1 ,x)f(T 2 ,x + Ax)d(Ax). (58) 



Inserting this into Eq. (|57| and taking into account that 
T\ is constant, and hence dT = dT 2 , we obtain 

f(s, x) = e- sV ^ At (f(s, x + Ax)) Ax , (59) 



s 



where the average is done over all realizations of Ax. 
With Taylor expansions around x of e ~ sV l x } At to the first 
order and of f(s, x + Ax) to the second order, Eq. (|59|) 
becomes 



f(s,x) = (l-sV[x]At) 

df(s,x) 



f(s,x) + 



Ox 



(Ax 



In a first order approach 



nD A 

Ax = — At- 
x 



l d 2 f(s,x) 
2 dx 2 



<2DAW t 



(Aa 



(60) 



(61) 



where AW t = Wt+At — W t , and thus, using Eq. ((2J, 

, . , nD . 

(Ax) = At. (62) 

x 

Then the mean value of (Ax 2 ) is, making again use of 
the zero mean property of the Wiener process, as well as 
of its autocorrelation function given in Eq. ([3]) , 



(Ax 2 ) = 2£>At + o(At). 



(63) 



Finally, putting V[X t ] = 1, we get the desired backward 
Kolmogorov equation for the first passage time PDF in 
Laplace space: 



d 2 f(s,x) ndf(s,x) 



dx 2 



x dx 



-f(s,x)=0. (64) 



C. Formulation of the boundary value problem 

We now proceed to the formulation of the boundary 
value problems corresponding to the solutions of the first- 
passage time PDFs, distinguishing between the three 
classes of boundaries the origin can belong to, as dis- 
cussed in Sec. IIIII On the right side we impose an ab- 
sorbing boundary at b: the first-passage time vanishes 
for x — > b~ , i.e., f{T,x — > b~) = S(T). Inserting this into 
Eq. ([M]) gives 



lim /(s,x) 

x— >b 



1. 



(65) 



The simplest case is if the zero level is an entrance 
boundary, i.e. n > 1. Starting from an inititial value 
Xq = x > 0, the zero level can never be reached, which 
corresponds to a reflecting wall at the origin. Applying 
standard arguments for reflecting boundaries [50J, the 
corresponding boundary condition is 



fen W«' a > = 0. 



>o+ 



dx 



(66) 



For n < — 1 the origin is an exit boundary. This means 
that it is impossible to reach any interior point of the 
state space if the initial point approaches the origin. This 



means that we have an absorbing boundary correspond- 
ing to 



lim f{s,x) 

x->0+ 



1, 



(67) 



and the first-passage time with respect to x = b will di- 
verge. Instead of the first-passage time the analysis of the 
previous section resulting in the backward Kolmogorov 
equation (l64l) together with the boundary conditions (|65|) 
and (j6"7| gives the first-exit time from the interval (0,6). 

In the case of a regular boundary, which happens for 
— 1 < n < 1, the behaviour is the most complicated. 
The process can both reach and leave the boundary zero, 
which means that also zero crossings are possible and the 
support of the process is the whole real axis. The first- 
exit time from (0, b) is again given by the same boundary 
condition problem as in the case of the exit boundary. 

For the sake of simplicity we rename /(s,x) =: y(x). 
Restricting the process to the positive half axis, our 
boundary value problem for the three kinds of bound- 
ary in the origin reads 



y"{x) 



-y'{x)- -pVix) = o, 

Ay(0) + By(o) = c, 



(68) 
(69) 



where 



^=1/(4)' B= lioJ- ^ 

An absorbing boundary at zero corresponds to 



1 




(71) 



whereas a reflecting boundary at zero corresponds to 

a = j ;; ;, i . c = [ v I . (72) 



Multiplying Eq. (J68|) with the integrating factor 
exp (J ^dx) leads to 



- (xV)' - -^ n y- (73) 

This is the canonical Sturm-Liouville form 1521 



(py'Y + mj = Awy, 



(74) 



with p(x) — x™, the weighting function w(x) = x", 
q[x) = 0, and the spectral parameter A = —s/D. 

We now observe that u := y — 1 transforms the ho- 
mogeneous problem (|68[) with inhomogeneous boundary 
conditions (|69p into an inhomogeneous problem with ho- 
mogeneous boundary conditions 



— (pu ) = Xwu + Xw, 
Au(0) + Bu(o) = 0, 



(75) 
(76) 



where u(a;) = (u(x), u'{x)) , and the two possible 
choices of A and c correspond to the Dirichlet problem 
and the Dirichlet-Neumann problem, respectively. This 
is easier to solve, since it determines a self-adjoint oper- 
ator C defined by 



Cu = — 
w 



■dm')' 



(77) 



in the weighted Hilbert space H = L 2 (J,w), where we 
have defined the open interval J = (0, b). This operator 
is not to be confused with the Laplace transformation 
operator Ct in Eq. (|54[) . which can be recognized from 
the index indicating the transformed variable. 

This can be seen as follows. Let u,v £ H; then the 
inner product is given by (u, v) = L uvw dx; taking into 
account that u and v satisfy the homogeneous boundary 
conditions (|76|) . we get after integrating twice by parts 



(u, Lv) = — I u{pv')' dx 
Jo 

r I 6 

= p(vu — uv ) 



o 



(pu')'v dx 



(Cu,v) 



(78) 



Using the definition from Eq. (|TT[) the boundary value 
problem given by Eqs. ([55HS5)) can be simplified to 



{C-Xl)u 
Au(0) + Bu(6) 



A, 
0. 



(79) 

(80) 



D. Formal solution of the boundary value problem 

We now exploit the property that the homogeneous 
boundary value problem with homogeneous boundary 
conditions 

{C-al)u = 0, (81) 

Au(0) + Bu(fc) = 0, (82) 

has nontrivial solutions Uk with eigenvalues a k , k G N, 



Cuk = otkUk- 



(83) 



Because C is self-adjoint, the eigenvalues a^ are real and 
the eigenfunctions Uk form an orthonormal basis of H. 
Furthermore a& > holds, since a k — (uk,Cuk). Hence 
the solution u of the inhomogeneous problem given by 
Eqs. ([79]) and (|80|) can be expressed through an expansion 
in this basis, 



y^CfcMfc, 
fe=i 



(84) 



with Ck — (uk , u) . Inserting u = 1 gives the normal- 
ization, X)feLi( u fci l) u fe = 1- The coefficients c& can be 
derived from Eq. d79l) : 

(85) 



{u kl Cu) - (u k ,Xu) = (u fe ,A). 



Again, making use of the definition of a self-adjoint oper- 
ator, we can pull C into the first component of the inner 
product. Employing Eq. (|83|) we get 



Ck 



(wfe,A) 



a.k - A 
The solution of the inhomogeneous problem reads 

v^ (wfe,A) 
u = 2_^ T u k- 



fc=i 



ak - A 



(86) 



(87) 



Because the eigenfunctions Uk do not depend on A = 
—s/D and the Laplace transformation is a linear opera- 
tion we obtain the inverse Laplace transform of y = 1 + u 
as 



y(T,x) = Cj^frxMT) 

oo 



fc=l 



oik - A 
5(T)+J2(ukA)u k [a k De- a » DT -S(T)] 



fc=i 



Since the Uk are normalized the two delta functions can- 
cel out. Returning to our original notation, we write the 
final result for the first-passage time (or first-exit time 
when appropriate) PDF as 



f{T 7 x) =Y,(uk,l)u k a k De- a ' 



DT 



(89) 



fc=i 



Of course, this PDF is normalized to 1: knowing that 
a k > we have 



/ f(T,x)dT = VVitfc,l)tt fc a:fcD / 
Jo k=l Jo 

OO 

= ^(u k ,l)uk = 1. 



,-ocuDT 



dT 



(90) 



fc=i 



We are now able to solve the specific boundary value 
problems for the three different kinds of boundaries at 
zero. 



E. Comparison of theory and simulation 

1. Simulation method 

To simulate the process X t that fulfills Eq. ([1]), we have 
used the Euler-Maruyama method [53 56], which in this 
case with an additive noise is identical to the higher- 
order Milstein method |54H56]. A generic autonomous 
stochastic differential equation 



dX t = fi(X t )dt + a(X t )dW t 



(91) 
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can be integrated between two successive times i„ and 
t n +i, giving 



X n +i — X n 



fjt(X t )dt + 



a(X t )dW t , (92) 



where X n is short for X tn . The approximation of the 
integrands to their value in t n , 






V(Xn), 



(93) 



yields the Euler method for the Ito case, called Euler- 
Maruyama [53j |. 



X n+1 = X n + fi(X n )At + a(X n )AW n , 



(94) 



where At = t n +i — t„ and AW n = W n+ i — W n ~ 
A(0, At) ~ v / AtiV(0, 1), i.e. AW n it is a normal random 
variable with PDF 



p(w) 



1 



V27rAt 



exp 



w 2 \ 
2At) 



(95) 



The Euler-Maruyama method has strong order of con- 
vergence 1/2. The Milstein method raises this to 1 
adding to the right-hand side of Eq. (|94| the correc- 
tion \a{X n )a'{X n )[{AW n f - At], where a'{X n ) = 
da(x)/dx\ x —x n - However, for an additive noise this 
derivative vanishes and so here the correction is zero. 
Schemes with order higher than 1 contain further terms 
some of which are nonzero also for additive noise, though 
many cancel out with respect to the general case, called 
multiplicative [44j . where a depends on X t . 
The approximation of the noise term as 



or as 



<?(X t ) 



°(Xt) 



a{X n ) + a(X n+1 ) 



X n + X n+ \ 



(96) 



(97) 



yields the corresponding method for the Stratonovich 
case; if a is continuous, both Eqs. (|96|) and (|97"1) lead 
to the same limit for At —¥ 0. This results in an implicit 
method, where to compute X n+ \ it is required to esti- 
mate it before; the predictor-corrector approach where 
X n+l in Eq. flHH) or flg7|) is approximated by Eq. (|M|) 
for the Ito case is known by the name of Euler-Heun or 
Heun 54, 55]. As already observed at the beginning of 
Sec. IIII Al with respect to the Fokker-Planck equation, 
both the Ito and the Stratonovich convention lead to the 
same result when the noise is additive as here. Inter- 
estingly the Milstein scheme represents both the order 1 
strong Ito- Taylor approximation and the order 1 strong 
Stratonovich- Taylor approximation, i.e. even in the mul- 
tiplicative case it coincides for both kinds of stochastic 
integral. 



In other words, the choice of X t within the discretiza- 
tion interval [X n ,X n+ i] affects the outcome of the inte- 
gration only as far as the dependence of the noise term 
a on X t is concerned, because the covariation of X t and 
of the Wiener process Wt driving the stochastic integral 
is not zero, [A t ,Wt] 7^ [57| (unfortunately closed in- 
tervals and covariations share the same notation). The 
choice of X t G [X ni X n+ i] has no influence on the in- 
tegration of the drift term \x with respect to t, and the 
choice of t G [t„,t n+1 ] does not matter for either fi or a 
if they depend on t, as [t,t] — and [t, W t ] = 0; in the 
three latter cases the same limit results for At — > 0. 

Eq. ([Ml) can be implemented straightforwardly in code. 
However, measuring the first-passage time with respect 
to a certain level needs a further refinement, since there is 
a finite hitting probability during each discretized time 
interval At, and thus the first-passage time is overesti- 
mated. An analytic expression for the probability that 
the process hits the level b during a discretization interval 
At was found by Mannella [58J . If we introduce the abre- 
viations [i n = fi(X n ), fib = fi{b) and fx' h = dfi{x)/dx\ x= b, 
the hitting probability reads 



P(hit) = exp I - 



A, 



2D (e 2 ^ At - 1) 



X n+1 -b+(X n - b)e^ 



At M6 

v 6 



-I 2 



1 



4L>At 



a-„ +i - [jr„ + Mn+ ^" +1 At 



(98) 



We can now summarize the simulation algorithm. We 
draw a Gaussian random number AW n using e.g. the 
Box-Muller method [59| and propagate the process X n 
by a time step At. If the propagated value exceeds the 
level b for the first time, i.e. X n +i > b, the process is 
terminated. Otherwise we check for missed hits in the 
discretization interval by drawing a uniformly distributed 
random number U G [0, 1) and accepting the hitting hy- 
pothesis if P(hit) > U; this fulfills the second terminating 
condition. In both cases the first-passage time is set to 
t n , i.e. the value before the propagation. 

In the case of an entrance boundary a further refine- 
ment of the simulation algorithm is possible. Knowing 
that the zero level can never be reached from the interior 
of the state space of the process, it is clear that negative 
values in the simulations must result from discretization 
errors. If this is the case we can reduce the time step 
until the propagated value of the process is positive. 



2. Entrance boundary 

As we know from the classification of the origin, we 
have an entrance boundary for n > 1 (i.e. v < 0). The 
general solution of the homogeneous differential equa- 
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tion (ED is 



t(x) = x v [AJ„ (y/^x) + BY V (y/ax)] 



(99) 



where J v and Y v are the Bessel functions of the first and 
second kind, respectively. 

Exploiting J' v (x) = J u -\ (x) — (y/x)J v (x) and an anal- 
ogous formula for Y v [46| , we obtain the derivative 



i'{x) — y/ax v [AJv_i 



BY V -! (Vax)] . (100) 



The relevant boundary conditions for f(s,x) given by 
Eqs. (|65p and (|66[) transform to u(b) — (absorption 
at X = b) and lim x _i.o u'{x) = (reflection at x — 0), 
respectively. 

To evaluate the eigenfunctions in the case of negative 
and integer v one can use the symmetry relation [46| 



j- v {z) = {-iyj v {z), 



(101) 



which holds for integer v 1 to see that the first term in 
Eq. (|100[) goes to zero for x — > 0, because its leading 
order term behaves as x. Since the Bessel functions of 
second kind diverge as x — > 0, the reflecting boundary 
condition can be fulfilled only if B = 0. 

The absorbing boundary condition at x = b determines 
the eigenvalues of the problem by the requirement that 
Jviyfoihb) = 0. Denoting the fcth zero of J v {x) by j k 
we thus have u k (x) = A k X v J v (jkx/b). The constant A}. 
is determined by the condition (u k ,Ui) = Ski- Remem- 
ber that the brackets denote the scalar product in the 
weighted Hilbert space with weigth w = x n . Observing 
the orthogonality relation (45| 



J„ 



) J, (jl ^)xdx = ^ b 2 J 2 u+1 (j k ) 8 M (102) 



Jk-r) -A 



one obtains Ak = v26 1 / Ju+i{jk), so that 

u k (x) = \f2b~ 1 x v J v (jk-r) /J v +i(jk)- 
We can further compute [43] 

V2b^ \ (jk/^y- 1 j v -x{j k 



(«fe,i) = 



3k 



T{v)J v+ i{j k ) J»+i(jk) 



(103) 



(104) 



For integer v a recurrence relation J u+ i(x) + J v —i(x) = 
2vJ v (x)/x holds, which, evaluated at the fcth zero of J„, 
delivers J v -i{jk) = —Jv+\{jk)- Hence Eq. (I104p simpli- 
fies to 



(uk, 1) = 



V2b l 



Jk 



tik/2) 



v-\ 



Y{v)J v+1 {j k ) 



(105) 



For non-integer v the Bessel functions J v and </_„ are 
two linear independent solutions of Eq. (|8ip . and it is 
more convenient to write the general solution as 



u(x) = x v \AJ V U/ax) + BJ- U [y/ax)\ . 



(106) 



Exploiting J'_ v (x) — —J\_ v (x)—(v/x)J- v {x) [4a], the 
derivative can be written as 

u'(x) = y/^x" [AJ v ^{y/ax) - BJi_„(y/ax)] . (107) 

In this case x i> J v -i(y/ax) diverges as x — > 0, and the 
left boundary condition requires A — 0. The right bound- 
ary condition determines the eigenvalues similarly as in 
the previous case; it is required that y/ab are the zeros 
j k of the Bessel function J_„. Again the second constant 
is evaluated using Eq. (J102I) . The normalized eigenfunc- 
tions are 



u k {x) = \/2 b l x u J- u [jk~\ /Ji- V (jk) 



with 



(w/c,l) 



V2b^j-\ 



(108) 



(109) 



For completeness we prove that the eigenfunctions 
given by Eqs. f)103|) and (|108[) fulfill the normalization 
condition (|9"0"|) . Setting x/b — z, for integer v Eqs. (|103[) 
and (fTUS]) yield 



El i\ »\T^ 2J v {j k z) 
(u k ,l)u k = z > — — - 
friJkJu+iijk) 

\2/j k f- v J v {j k z) 



fc=l 



{Vjk) l - V 






T(v)J v+ i(j k ) 

2J v {j k z) 



T { v )Ju+\{Jk) jkJv+l{ik) 



1, (no) 



where we have used the Fourier-Bessel expansions 
v ^ 2J v (j k z) 



f-[jkJv+l{jk) 



(2/j k ) 2 -»J u ( Jk z) 



2,1 z ' ~ Z^ r,,,\ r-> 



fc=i 



mJi+iUk) 



(in) 



(112) 



Eq. (Illip is found in Watson 45] , Eq. (|112[) is proved in 
the appendix. For non- integer v Eqs. (|108p . (|109p and 
(fTTTj) yield 



2J-„(j k z) 



V(u fe ,l)u fe =z v 

f x f^jkJi-Ajk) 



(113) 



Fig. Q] shows the analytical results obtained with Wol- 
fram Mathematica 7.0 by truncating the sum in Eq. (|59")l 
after the first 200 terms, and normalized histograms gen- 
erated with 10 million simulation runs done as explained 
in Sec. HVETl The agreement is perfect. The CPU time 
needed for an analytical curve is a few seconds, while 
that for a histogram with 10 million runs, which is the 
number used for all histograms in this paper, ranges from 
a few minutes to two days, depending on the time step, 
the starting position, and the upper boundary b. We 
used the rani uniform random number generator [601 ] 
and the GNU C++ compiler (g++) version 4.1.2 with 
the -03 optimization option on a 2.2 GHz AMD Athlon 
64 "Winchester" processor with Fedora Core 7 Linux. 
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n = 2.5, x = l,b = 4 
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FIG. 1: First passage time PDF f(Tt) when the origin is an 
entrance boundary. The parameter n, the starting position x, 
and the upper boundary b are given in the insets; the diffusion 
coefficient D is 1. The analytical results (lines) are perfectly 
covered by the normalized histograms obtained from simula- 
tion (circles). 



FIG. 2: First exit time PDF f(T 0i b) when the origin is an 
exit boundary. The parameter n, the starting position x, and 
the upper boundary b are given in the insets; the diffusion 
coefficient D is 1. The analytical results (lines) are perfectly 
covered by the normalized histograms obtained from simula- 
tion (circles). 



3. Exit boundary 

For n < — 1, i.e. v > 1, the zero level is an exit bound- 
ary and it is impossible to reach any interior point b, pro- 
vided that the starting point of the process is sufficiently 
close to the boundary. This, repeating the arguments 
of Sec. IIIIA1 corresponds to a collapse of the PDF to a 
delta function 5{x) in a finite time. Hence in general the 
first-passage time with respect to x — b > will diverge. 

However, with the absorbing boundary at the origin, 
where naturally w(0) = 0, and imposing an absorbing 
boundary condition at the upper boundary, u(b) = 0, we 
have a boundary value problem with a solution that is 
the PDF of the first-exit time To.b = min{Xb,Tf,} from 
the interval (0,6). 

The first terms in Eqs. (|M| or (|106[) vanish in the limit 
x — > 0, whereas the second terms diverge. Therefore the 
constant B must be zero in order to fulfill u{b) = 0. 
As in the case of an entrance boundary, the eigenvalues 
are determined by the condition J v (yfctkb) — and the 
eigenfunctions are given by Eq. (1103J) . 

In Fig. [2] the theoretical curves are again compared 
with the results obtained numerically. It is interesting to 
notice that the PDF is bimodal for a range of parameters. 



by Eq. (|103[) . which was also the result for the entrance 
boundary in the case of negative and integer v. Fig. [3] 
shows that the theoretical results agree perfectly with the 
simulations. 



n;„„) 



n = -0.5, x = 4,b = 5 




FIG. 3: First exit time PDF /(Tb,i>) when the origin is a regu- 
lar boundary imposed to be absorbing. The parameter n, the 
starting position x, and the upper boundary b are given in the 
insets; the diffusion coefficient D is 1. The analytical results 
(lines) are perfectly covered by the normalized histograms ob- 
tained from simulation (circles). 



4- Regular boundary 

For — 1 < n < 1, i.e. < v < 1, the origin is a regular 
boundary, and in accordance with Karlin and Taylor [48| 
it is possible to impose different boundary conditions in 
a consistent way. 

Imposing an absorbing boundary condition at the ori- 
gin gives the PDF of the first-exit time. The eigenval- 
ues are computed in the same way as for an exit bound- 
ary at the origin, and the eigenfunctions are again given 



It is also possible to impose a reflecting boundary con- 
dition at the origin, and hence the eigenfunctions are 
computed in the same way as in the case of an entrance 
boundary, of course inserting the respective value of n. 
For a few values of n in the range < n < 1 we have 
compared the theoretical curves resulting from the as- 
sumption that the origin is reflecting with histograms 
from simulations where we have allowed zero crossings; 
see the squares in Fig. @] It appears that for small times 
there is a good agreement, whereas for larger times there 
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are differences: the maxima of the histograms are higher 
than predicted by the theory, and the tails obtained by 
simulations are flatter than the theoretical tails. So we 
can clearly conclude that the origin is not naturally re- 
flecting for this range of n, but, as one can see in Fig. |4l 
total reflection is approached when n is approaching the 
limit where the origin is an entrance boundary, namely 

71 = 1. 



P(x,3) 



n = 0.8, x= l,b = 4 




n = 0.5 



• At = 10" 1 




*•* 


■ At = 1(T 2 




1 1 


At = 1(T J 




> 1 


• At = 1(T 4 






At = 1(T 5 


0.1 


- : 


-' 







^.mni?" 



p(x,3) 



n= 1.0 

• At= 10 _1 

■ At = 1(T 2 
At = 10~ 3 

■ At = 1(T 4 
At = 10~ 5 



.» 



FIG. 4: First passage time PDF f(Tb) when the origin is 
a regular boundary imposed to be reflecting. The parame- 
ter n, the starting position x, and the upper boundary b are 
given in the insets; the diffusion coefficient D is 1. The ana- 
lytical results (lines) are perfectly covered by the normalized 
histograms obtained from simulation (circles). Mismatching 
simulation results (squares) arise if zero crossings are allowed. 

To explain this phenomenon we recall what we have 
mentioned earlier: one might think intuitively that zero 
crossings are not possible for non-zero values of n since 
the drift term explodes near the origin, and the latter 
either reflects or absorbs the process for all times. How- 
ever, applying Feller's formal classification scheme one 
can see that zero-crossings are actually allowed for a reg- 
ular boundary at the origin, i.e. — 1 < n < 1. This is 
further confirmed by Fig. [5] 

Knowing this and the fact that according to Karlin 
and Taylor [48j a process can spend a finite time in the 
vicinity of a regular boundary, we can explain the plots 
qualitatively. The paths that are able to escape the in- 
fluence of the origin will quickly hit the boundary b fol- 
lowing the same rules as for the entrance boundary; they 
are basically driven by the drift term. The deviation in 
the tail of the PDFs is due to the positive amount of 
time spent in the vicinity of the origin, which is called 
the sticky boundary phenomenon [481 ] . and to the mul- 
tiple zero crossings. Fig. [6] shows a logarithmic plot of 
the first-passage time PDFs obtained by simulations, and 
one can see that the latter are heavy-tailed, i.e. they ex- 
hibit a power-law decay for long times. This is in contrast 
to the exponential decay obtained for the other types of 
boundary. 

However, there is a good agreement between theory 
and simulation if we impose a reflecting origin in the 
simulation too, meaning that we consider the origin as a 
hard reflecting wall; see the circles in Fig. @] 

For — 1 < n < the first-passage times diverge if we 
do not impose any artificial boundary condition, since 
the drift term is always negative if X t > and positive if 
X t < 0, meaning that the process is always attracted, but 
not totally absorbed, by the origin. On the other hand, 






FIG. 5: PDFs of the stochastic process X(t) at t — 3 from 
simulation with different integration time step At for n = 0.5 
(left) and n — 1.0 (right); the starting position xq and the 
diffusion coefficient D are 1. For n — 0.5 the peak in the 
negative domain increases with decreasing time step, whereas 
for n — 1 it decreases. This suggests that in the latter case, 
where the origin is an entrance boundary, the zero crossings 
are an artifact due to the discretization of time, whereas in 
the former case, where the origin is regular, the zero crossings 
are genuine. 



f(T b ) 

ir 







n = 0.5, x= l,b = 4 
n = 0.8, x= l,b = 4 



50 100 201 



ifl< 



FIG. 6: Normalized histograms of first-passage times with 
respect to the level b obtained by simulation when the origin is 
a regular boundary. The parameter n, the starting position x, 
and the upper boundary b are given in the inset; the diffusion 
coefficient D is 1. The tails can be fitted by power laws with 
exponents -1.11 and -1.17, respectively. 



imposing total absorption at the origin corresponds to 
the computation of the PDF of the first-exit times as 
shown in Fig. [3] 

It is interesting to note that the case n = 0, i.e. the 
Wiener process, belongs to this class. The origin is not 
a singular point and the first-passage time PDF with re- 
spect to a level b starting at xq is given by Eq. Q , which 
for long times is a power law with exponent —3/2. 
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CONCLUSIONS 



resulting in 



We have computed first-passage and first-exit time 
PDFs for a stochastic process with applications in many 
physical, chemical, biological, economical and financial 
problems. Depending on the nature of the boundary 
at the origin, we have found analytical solutions for the 
first-passage and first-exit time PDFs for all cases, ex- 
cept for the first-passage time PDF in the case of a reg- 
ular boundary at the origin. In the latter case we have 
found an analytical solution for the first-exit time PDF 
and approximations for the first-passage time PDF for 
short times. For this specific stochastic process regularity 
of the boundary at zero can include behaviours ranging 
from total absorption to total reflection, with intermedi- 
ate behaviours like elastic and sticky boundaries [48j . It 
is interesting that sticky boundaries may be applied e.g. 
to simulate the partial adsorption of polymer molecules 
to walls and for the modeling of solvent quality [6l| . In 
possible future projects this could be investigated more 
thoroughly and regarded from the perspective of interac- 
tions between molecules and boundary surfaces, which is 
closely connected to another project of two of us [62, [63| , 
where discotic liquid crystals confined in cylindric geome- 
tries [64J are studied via molecular dynamics simulations. 



Ck 



Ju+lUk) Jo 

For f(z) = Z~ v 



l 21 

J v {jkz)f{z)zdz = -2 — — -T-. (116) 

J v+lUk) 



Ik = / Ju(jkz)z " ' dz. 
Jo 



In order to exploit the equation |6£ 



(117) 



J v {z)z l ~ v dz = -J v ^{z)z x - v , (118) 



we substitute jkz — a and get 



:i k 



" J ' J v {oi)a 1 - v da 



Jk 2 [~Ju-i(a)a 



1— 1/1 Ok 





Jk- 2 



lim J v _ 1 {a)a 1 ~ v - J v -i{jk)Jk (- 119 ) 
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Appendix 

We prove the Fourier-Bessel expansion given in 
Eq. (|112[) . If a function f(z) is represented in an orthog- 
onal basis of Bessel functions of the first kind J v (jkz), 
where jk is the fcth zero of J u {z), i.e. J v (jk) = 0, 



f{z) = ^c k Ju{Jkz), 



(114) 



The limit is 



lim 



(-1)' /«N 



ima'^i^fa) = lim a 1 "" V „ \, ' : (-) 



21+v-l 



E 

1=0 
2 l~u 



T(v) 



(-1)'2 



l nl-v-2l 



lim a 



2/ 



l\T{l + v) a^O 



yielding 



I k=.] k 



Ju-l(jk) 



r M fk 



(120) 



(121) 



fc=i 



and the orthogonality relation is given by Eq. (|102j) . the 
Zth coefficient q can be obtained from the scalar product 
of f{z) with the Ith basis set element J u (jiz), 

/ J v {Jiz)f{z)zdz = y^Cfc / J v (j k z)J v (jiz)zdz 

Jo ,. , Jo 



^4-Jl+iiik)^ 



II, 



fe=l 



■£-%+i(Ji)> 



(115) 



Thus 



E 
fc=i 



(2/jfc) 



2J„-l(j/c) 



T{v) 



Jk 



JyjjkZ) 
Jt+ltikY 



(122) 



Subtracting the Fourier-Bessel expansion of z v ', 
Eq. (jllll) . the second term in square braces 
cancels out because of the recurrence identity 
Ju-i(jk) + J v +xijk) = 2vJ v {jk )/jk = that we 
have already used to simplify Eq. (|104p to Eq. (|105[) . 
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